pcvrDDPSC Data Science Core, August 2023
Josh Sumner
pcvr GoalsgrowthSSfitGrowthgrowthPlot/Model Visualizationbrms directlypcvr GoalsCurrently pcvr aims to:
There is room for goals to evolve based on feedback and scientific needs.
Pre-work was to install R, Rstudio, and pcvr with dependencies.
plantCV allows for user friendly high throughput image based phenotyping
Resulting data follows individuals over time, which changes our statistical needs.
Longitudinal Data is:
Bayesian modeling allows us to account for all these problems via a more flexible interface than frequentist methods.
Bayesian modeling also allows for non-linear, probability driven hypothesis testing.
In a Bayesian context we flip “random” and “fixed” elements.
| Fixed | Random | Interpretation | |
|---|---|---|---|
| Frequentist | True Effect | Data | If the True Effect is 0 then there is an \(\alpha\cdot100\)% chance of estimating an effect of this size or more. |
| Bayesian | Data | True Effect | Given the estimated effect from our data there is a P probability of the True Effect being a difference of at least X |
Before moving on we should note that each group in your data will have parameters fit to it. If you have many groups then fitting models to only a few groups at a time is likely to make your life easier. Larger models will fit, but they can take a very long time.
It’s generally been easier for the authors to fit 50 models each with 4 groups in them then compare across models than to fit 1 model with 200 groups in it.
There are 6 main growth models supported in pcvr.
3 are asymptotic, 3 are non-asympototic.
There are an additional 3 sigmoidal models based on the Extreme Value Distribution. Those are Weibull, Frechet, and Gumbel. The authors generally prefer Gompertz to these options but for your data it is possible that these could be a better fit.
There are also two double sigmoid curves intended for use with recovery experiments.
Generalized Additive Models (GAMs) are supported but discouraged due to poor interpretability.
Intercept only models are supported, although they are only intended for use in segmented models or to represent homoskedasticity as a sub-model.
Segmented growth models are also supported using “model1 + model2” syntax in both growthSim and growthSS.
For now we will focus on single models but there will be examples of segmented models later on.
Finally, the “decay” keyword can be used to specify a decay model instead. For now we will focus on growth models as these constitute the large majority of models used in plant phenotyping.
Survival models can also be specified using the “survival” keyword. Using the brms backend these models use the weibull distribution by default but “binomial” family can also be specified by model = "survival binomial".
For details please see the growthSS documentation.
growthSSAny of the aforementioned models can be specified easily using growthSS.
growthSS - formWith a model specified a rough formula is required to parse your data to fit the model.
The layout of that formula is:
outcome ~ time|individual/group
growthSS - form 2Here we would use y~time|id/group
growthSS - form 3Generally it makes sense to visually check that your formula covers your experimental design.
Note that it is fine for id to be duplicated between groups, but not within groups
growthSS - sigmaRecall the heteroskedasticity problem, shown again with our simulated data:
growthSS - sigmaThere are lots of ways to model a trend like that we see for sigma.
pcvr allows any of the model options in growthSS to also be applied to model variance.
growthSS - Intercept sigma| Pros | Cons |
|---|---|
| Faster model fitting | Very inaccurate intervals at early timepoints |
growthSS - Linear sigma| Pros | Cons |
|---|---|
| Models still fit quickly | Variance tends to increase non-linearly |
| Easy testing on variance model |
growthSS - Linear sigmagrowthSS - Gompertz sigma| Pros | Cons |
|---|---|
| Models fit much faster than splines | Slightly slower than linear sub-models |
| Variance is often asymptotic | Requires priors on sigma model |
| Easy testing on variance model |
Note that these traits are broadly true of logistic and monomolecular sub models as well.
growthSS - Gompertz sigmagrowthSS - Spline sigma| Pros | Cons |
|---|---|
| Very flexible and accurate model for sigma | Significantly slower than other options |
| Fewer priors | Splines can be a black-box |
growthSS - Spline sigmagrowthSS - other sigma modelsYou can always add a new sigma formula if something else fits your needs better.
growthSS - Distributional ModelsHere we have limited the examples to talk about sigma, a parameter of the Student T distribution that our model belongs to. Written another way we might say all the previous methods are modeling:
Y ~ T(mu ~ main growth formula, sigma ~ sigma formula, nu ~ 1)
growthSS - Distributional Models 2In pcvr the Student T family is the default for these models, but other distributions are supported through "distribution: model" syntax.
In general this is only for special cases where the Gaussian/T does not capture some important quality of the data. An obvious example could be leaf counts, which might be modelled as "poisson: monomolecular", for instance.
growthSS - Distributional Models 3You can specify sigma as a list of formulas to model different distributional parameters separately. Most of the time this is overkill and adds unnecessary complexity, but the option exists for certain cases such as ZINB where modeling mean, shape, and zero inflation per group may make sense.
For details on supported families and their parameterization try running growthSS(..., sigma=NULL,...) and examining the priors, or checking ?brmfamily.
growthSS - priorsBayesian statistics combine prior distributions and collected data to form a posterior distribution.
Luckily, in the growth model context it is pretty easy to set “good priors”.
growthSS - priors“Good priors” are generally mildly informative, but not very strong.
They provide some well vetted evidence, but do not overpower the data.
growthSS - priorsFor our setting we know growth is positive and we should have basic impressions of what sizes are possible.
At the “weakest” side of these priors we at least know growth is positive and the camera only can measure some finite space.
growthSS - priors 2Default priors in growthSS are log-normal
\(\text{log}~N(\mu, 0.25)\)
This has the benefit of giving a long right tail and strictly positive values while only requiring us to provide \(\mu\).
growthSS - priors 3We can see what those log-normal distributions look like with plotPrior.
growthSS - priors 3growthSS - priors 4Those distributions can still be somewhat abstract, so we can simulate draws from the priors and see what those values yield in our growth model.
growthSS - priors 4growthSS - priors 5Our final call to growthSS will look like this for our sample data.
fitGrowthNow that we have the components for our model from growthSS we can fit the model with fitGrowth.
fitGrowth 2This will call Stan outside of R to run Markov Chain Monte Carlo (MCMC) to get draws from the posterior distributions. We can control how Stan runs with additional arguments to fitGrowth, although the only required argument is the output from growthSS.
fitGrowth 2Here we specify our ss argument to be the output from growthSS and tell the model to use 4 cores so that the chains run entirely in parallel, but the rest of this model is using defaults.
fitGrowth 2Note that there are lots of arguments that can be passed to brms::brm via fitGrowth.
One that can be very helpful for fitting complex models is the control argument, where we can control the sampler’s behavior.
fitGrowth 2adapt_delta and tree_depth are both used to reduce the number of “divergent transitions” which are times that the sampler has some departure from the True path and which can compromise the results.
fitGrowth 3fitGrowth returns a brmsfit object, see ?brmsfit and methods(class="brmsfit") for general information.
Within pcvr there are several functions for visualizing these objects.
growthPlotgrowthPlot can be used to plot credible intervals of your model.
These plots can show one of the benefits of an asymptotic sub model well.
Here we check our model predictions to 35 days.
And now we check those predictions from a spline model, where the basis functions are not suited for data past day 25.
We can also plot the posterior distributions and test hypotheses with brmViolin.
Here hypotheses are tested with brms::hypothesis.
brms::hypothesis allows for incredibly flexible hypothesis testing.
Here we test for an asymptote for group A at least 20% larger than that of group B.
Segmented Models are specified using “model1 + model2” syntax, with “+” representing a change point.
Currently only two phases (one changepoint) are recommended. While more will work they will slow down significantly and may require more fine tuning.
These segmented models can also be used to specify sub-models of sigma.
simdf<-growthSim(model = "linear + linear",
n=20, t=25,
params = list("linear1A"=c(15, 12), "changePoint1"=c(8, 6), "linear2A"=c(3, 5) ))
ss<-growthSS(model = "linear + linear", form=y~time|id/group, sigma="spline",
start = list("linear1A"=10, "changePoint1"=5, "linear2A"=2 ),
df=simdf, type = "brms")
fit <- fitGrowth(ss, backend="cmdstanr", iter=500, chains=1, cores=1)
growthPlot(fit=fit, form=ss$pcvrForm, df = ss$df)simdf<-growthSim("linear + logistic", n=20, t=25,
params = list("linear1A"=c(15, 12), "changePoint1"=c(8, 6),
"logistic2A"=c(100, 150), "logistic2B"=c(10, 8),
"logistic2C"=c(3, 2.5) ))
ss<-growthSS(model = "linear + logistic", form=y~time|id/group, sigma="spline",
list("linear1A"=10, "changePoint1"=5,
"logistic2A"=100, "logistic2B"=10, "logistic2C"=3),
df=simdf, type = "brms")
fit <- fitGrowth(ss, backend="cmdstanr", iter=500, chains=1, cores=1)
growthPlot(fit=fit, form=ss$pcvrForm, df = ss$df)simdf<-growthSim("linear + linear + linear", n=25, t=50,
params = list("linear1A"=c(10, 12), "changePoint1"=c(8, 6),
"linear2A"=c(1, 2), "changePoint2"=c(25, 30),"linear3A"=c(20, 24)))
ss<-growthSS(model = "linear + linear + linear", form=y~time|id/group, sigma="spline",
list("linear1A"=10, "changePoint1"=5,
"linear2A"=2, "changePoint2"=15,
"linear3A"=5 ), df=simdf, type = "brms")
fit <- fitGrowth(ss, backend="cmdstanr", iter=500, chains=1, cores=1)
plot <- growthPlot(fit=fit, form=ss$pcvrForm, df = ss$df)ss<-growthSS(model = "int + int", form=y~time|id/group, sigma="int + int",
list("int1"=10, "changePoint1"=10, "int2"=20, #main model
"sigmaint1"=10, "sigmachangePoint1"=10, "sigmaint2"=10), # sub model
df=simdf, type = "brms")
fit <- fitGrowth(ss, backend="cmdstanr", iter=500, chains=1, cores=1)
plot <- growthPlot(fit=fit, form=ss$pcvrForm, df = ss$df)ss<-growthSS(model = "int + linear", form=y~time|id/group, sigma="int + linear",
list("int1"=10, "changePoint1"=10, "linear2A"=20,
"sigmaint1"=10, "sigmachangePoint1"=10, "sigmalinear2A"=10),
df=simdf, type = "brms")
fit <- fitGrowth(ss, backend="cmdstanr", iter=500, chains=1, cores=1)
plot <- growthPlot(fit=fit, form=ss$pcvrForm, df = ss$df, timeRange = 1:40)ss<-growthSS(model = "int+logistic", form=y~time|id/group, sigma="int + spline",
list("int1" = 5, "changePoint1"=10 ,
'logistic2A' = 130, 'logistic2B' = 10, "logistic2C" = 3,
'sigmaint1' = 5, "sigmachangePoint1"=15),
df=simdf, type = "brms")
fit <- fitGrowth(ss, backend="cmdstanr", iter=500, chains=1, cores=1)
plot <- growthPlot(fit=fit, form=ss$pcvrForm, df = ss$df)df <- growthSim("logistic", n=20, t=25,
params = list("A"=c(200,160), "B"=c(13, 11), "C"=c(3, 3.5)))
ss <- growthSS(model = "survival weibull", type="brms",
form = y > 100 ~ time|id/group,
df = df, start=c(0,5))
fit <- fitGrowth(ss,iter = 600, cores = 2, chains = 2, backend = "cmdstanr")
plot <- growthPlot(fit=fit, form=ss$pcvrForm, df = ss$df)Here the input data is standard phenotype data with a cutoff to represent the event (germination for instance) on the left hand side of the formula.
df <-growthSim("count: logistic", n=20, t=25,
params = list("A"=c(10,12), "B"=c(13, 11), "C"=c(3, 3.5)))
ss <- growthSS(model = "poisson: logistic", # specify poisson family
form = y ~ time|id/group,
sigma=NULL, # poisson only has one parameter
df = df, start=list("A"=8, "B"=10, "C"=3) )
fit <- fitGrowth(ss,iter = 2000, cores = 4, chains = 4)
plot <- growthPlot(fit=fit, form=ss$pcvrForm, df = ss$df)
ggsave("~/pcvr/tutorials/pcvrTutorial_advancedGrowthModeling/exCountModel.png", plot, width=6,height=4, dpi=300, bg="#ffffff")John Kruschke wrote a paper on the Bayesian Analysis and Reporting Guidelines (BARG) to aid in transparency and reproducibility when using Bayesian methods.
In pcvr some of what Kruschke recommends can be accessed from a model/list of models using the barg function, see ?barg for details on what that entails.
brms directlyThese functions are all to help use common growth models more easily.
The choices in pcvr are a small subset of what is possible with brms, which itself is more limited than Stan.
brms directlyOur gompertz sigma model looks like this in brms:
prior1 <- prior(gamma(2,0.1), class="nu", lb=0.001)+
prior(lognormal(log(130), .25),nlpar = "A", lb = 0) +
prior(lognormal(log(12), .25), nlpar = "B", lb = 0) +
prior(lognormal(log(1.2), .25), nlpar = "C", lb = 0)+
prior(lognormal(log(25), .25),nlpar = "subA", lb = 0) +
prior(lognormal(log(20), .25), nlpar = "subB", lb = 0) +
prior(lognormal(log(1.2), .25), nlpar = "subC", lb = 0)
form_b <- bf(y ~ A*exp(-B*exp(-C*time)),
nlf(sigma ~ subA*exp(-subB*exp(-subC*time))),
A+B+C+subA+subB+subC ~ 0+group,
autocor = ~arma(~time|sample:group,1,1),
nl = TRUE )
fit_g2 <- brm(form_b, family = student, prior = prior1, data = simdf,
iter = 1000, cores = 4, chains = 4, backend = "cmdstanr", silent = 0,
control = list(adapt_delta = 0.999,max_treedepth = 20),
init = 0 ) # chain initialization at 0 for simplicitybrms directlyIt can be more work to try new options in brms or Stan, but if you have a situation not well represented by the existing models then it may be necessary.
If you run into a novel situation please reach out and we will try to come up with a solution and add it to pcvr if possible.
Good ways to reach out are the help-datascience slack channel and pcvr github repository.